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We introduce a method of molecular dynamics in shape 
space aimed at metal clusters. The ionic degrees of freedom 
are described via a dynamically deformable jellium with in- 
ertia parameters derived from an incompressible, irrotational 
flow. The shell correction method is used to calculate the elec- 
tronic potential energy surface underlying the dynamics. Our 
flnite-temperature simulations of Agi4 and its ions, follow- 
ing the negative-to-neutral-to-positive scheme, demonstrate 
the potential of pump-and-probe ultrashort laser pulses as a 
spectroscopy of cluster-shape vibrations. 



The study of the motion of an incompressible and ir- 
rotational fluid mass with an ellipsoidal surface has its 
origins in investigations on the theory of the shape of the 
Earth. Such studies started with Newton and Maclaurin 
and were continued by others, in particular for the case 
of a varying surface by Dirichlet, Dedekind, and Rie- 
mann (see Ref. |]). In these early studies in the realm of 
astronomy and astrophysics, the forces were due to the 
gravitational field. More recently, similar ideas were used 
in nuclear physics in. investigations of low-energy collec- 
tive isoscalar modesH (as well as of the fission processB in 
atomic nuclei). Additionally, many feai-ures of the anal- 
ysis of these nuclear collective modeso go back to the 
work of Rayleigh on the normal modes of a classical liq- 
uid drop.Q 

In this letter, we introduce and apply a method of 
molecular dynamics in shape space (MDSS) aimed at 
studies of materials clusters, with the potential energy 
surfaces (PES's), which underlie the dynamics, deter- 
mined via electronic struct w-e, calculations using the shell 
correction method (SCM)Jj~Q and with inertia parame- 
ters specified according to the aforementioned theory of 
an irrotational fluid mass. 

In the past several years, first-principles molecular- 
dynamics (FPMD) simulation methods, where the dy- 
namics of the ions is evaluated on concurrently calculated 
electronic PES's [e.g., through local density functional 
(LDP) theory], have been introduced and emnLayed in 
investigations of atomic and molecular clusters B-EEl The 
MDSS method complements such techniques. Rather 
than simulating the dynamical trajectories of the indi- 
vidual ions as in other FPMD methods, the MDSS fo- 
cuses opithe dynamics of collective shape variations of 
clustersO (with the ionic degrees of freedom described 
via a dynamically deformable jellium), allowing for reli- 



able and economical simulations in a broad cluster size 
range. 

Over the past several years, the development and ap- 
plication of fast spectroscopical techniques opened new 
avenues for probing dynamical processes in matter (both 
in the molecular and condensed-phase regimes) with 
unprecedented temporal resolution (in the femtosecond 
rangell3). As an illustration, we apply the MDSS to an 
investigation of the shape dynamics of a small neutral 
cluster (Agi4), motivated by recent experiments using 
pump-and-probe ultrashort lases, pulses. In these exper- 
iments, referred to as NeNePoJl3 negative-to-neutral-to- 
positive, mass selected neutral clusters are prepared by 
electron vertical detachment from the associated nega- 
tive ions, and their subsequent time evolution is probed 
by photoionization, which can be time-delayed with ref- 
erence to the initial detachment. Our finite-temperature 
MDSS simulations of Agi4 and of associated ions, follow- 
ing the NeNePo sequence, demonstrate the potential of 
such experiments as a spectroscopy of cluster shape vi- 
brations. As the jellium description of metal clusters is 
known to be most adequate above the Debye tempera- 
ture, 6_D, we demonstrate the MDSS method for cluster 
shape dynamics at T > 9£)(Ag) (= 215 K). 

Modelling the dynamics of the ellipsoidal jellium back- 
ground as that of an irrotational, incompressible fluid, its 
kinetic energy is givenlll by 

r(a, h) = Avrpa'fe'c'(a' ^+b'^+c'^) , (1) 
to 

where p is the mass density of the fluid, and a', and c' 
are the principal semi-axes of the ellipsoid. Since the PES 
V{a' , b') depends only on two semi-axes, we can eliminate 
the third semi-axis c' from the kinetic energy by using 
the volume conservation, namely a'b'c' = R^. By further 
defining reduced semi-axes a = a' / R, b = b' / R, and c = 
c'/i?, we can rewrite the kinetic energy as, 

T (a, b) = ^Maaa^ + Mabab + ^Ahbb^ , (2) 

where the elements of the inertial-mass matrix M are 
Maa = Cil + 1/a^b^), Mab = A/ft, = C / a^b\ and Ahb - 
C(l 4- lla%^), with C = 47ri?V/15. 

At constant energy, the motion of the shape will be 
governed by Hamilton's equations associated with the 
hamiltonian H — T(a, b)+V{a, b). Introducing new sym- 
bols a ^ qi, b ^ q2, Pa ^ Pi, and pb ^ P2, for the 
generalized coordinates a and b and for the associated 
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canonical momenta pa and pb, Hamilton's equations of 
motion are given by 

2 

g.=^(A/-i),,p, , (3) 

where M^^ is the inverse of the inertial-mass matrix 
M . These equations are solved using a 6th-order Runge- 
Kutta method. 

MDSS simulations can be carried out at constant en- 
ergy or constant temperature. The latter (canonical- 
ensemble) simulations can be achieved through a 
straightforward modification of Hamilton's equations ajG-, 
cording to the Nose-Hoover-thermostat prescription,til 
namely by adding a variable friction term —PiTj/Q to the 
right-hand-side of the pi equation, where Q is a constant 
characterizing the reservoir and f] = 2T{qi^q2) — 2k^T, 
with T the temperature and k-Q the Boltzmann constant. 
Since, however, the triaxial MDSS-involves two indepen- 
dent variables only, we use a chainEj of three Nose- Hoover 
thermostats in order to guarantee that the equilibrium 
ensemble will exhibit canonically distributed positions 
and momenta. 

The SCM approach has been derived by us originally 
from the LDF theory (for details, see Re|^||,0). This 
method, as well as a semiempirical versionErLI of it which 
we employ here to determine the PES's [herir^ ®\^ 
has been shown to yield accurate resultsQ^la'Ej when com- 
pared to experiments and/or self-consistent Kohn-Sham 
LDF calculations (when available) for several metal- 
cluster systems. Particularly pertinent to this study are 
SCM results for the vertical Ionization Potential (vIP) of 
Agi4 and the electron vertical Detachment Energy (vDE) 
of Agj"4; the SCM results, calculated for the equihb- 
rium configurations, are vIP=5.84 eV and vDE=2.20 eV, 
in jescellent agreement with measured values, vIP=5.90 
eVllZl and vDE ~ 2.2 eV.El 

We prepare first an equilibrium ensemble of 10'^ phase- 
space points [(gi)o, (92)0, (pi)o, and (^2)0] of the an- 
ion (Agj"^) by selecting them from a long trajectory (up 
to 400 ps) on the PES of Agj"4 generated using the 
Nose-Hoover dynamics at a given temperature T . These 
phase-space points are subsequently used as initial val- 
ues for geneLating shorter (up to 3 ps), constant-euergy 
trajectories,!!^ with an integration time step of 0.5 fs,E2l on 
the PES of the neutral Agi4 (which is produced from the 
anion by electron photodetachment caused by the first 
laser pulse). At each time step, we determine the clas- 
sical probability, P(i), that vertical ionization of the dy- 
namically evolving neutral cluster is energetically allowed 
through the two-photon absorption process induced by 
the probe-laser pulse. V{t) is that fraction of the in- 
stantaneous shape configurations (averaged over the 10'^ 
constant-energy trajectories) with a vertical ionization 



potential vIP< 2Tujj, where huj is the single-photon en- 
ergy. 

The probability Vit) is a classical quantity which in 
principle is not the sole factor determining the yield of 
the photoionization process (e.g., the photoion current of 
the cationic daughters). Indeed the photoion current is 
proportional to the product of the classical V{t) with the 
quantum mechanical, vibrational, Frank-Condon overlap 
factor .Ell If the neutral clusters happen to be produced in 
shape configurations far away from the PES minima of 
the positive ion, the Franck-Condon overlap factor can 
be so low initially that virtually no positive ions are gen- 
erated, regardless of the value of V{t), and thus the onset 
of the photoion current can exhibit a characteristic delay. 
The NeNePft-fBtperiment for the Silver trimer conforms 
to this caseJl3'E2l However, the case of Agi4 is different, 
since here the PES minima of the neutral and the positive 
clusters overlap substantially (see Fig. 1). As a result, in 
the latter case, the Franck-Condon factor for the ioniza- 
tion process is close to unity, and we can safely assume 
that the photoion current is mainly proportional to the 
classical probability V{t). 

In Fig. 1, we display the PES's for Agj^^ (bottom 
panel), Agi4 (middle), and Ag^4 (top). For the neutral 
Agi4, which is axially symmetric, there are two energet- 
ically distinct isomers - an oblate (O) one and a prolate 
(P) one. On the associated PES (Fig. 1, middle panel), 
there are six local minima, three (energetically degener- 
ate ones) corresponding to the oblate 0-isomer (having 
two equal axes larger than unity; the case a = b > 1 
with c = \/a? is easily identifiable; the other two follow 
by circular permutation of the indices a, 6, and c, i.e., 
pseudorotations) and the other three (energetically de- 
generate ones) to the P-isomer (having two equal axes 
smaller than unity; the case of a = 6 < 1 with c = 1 /a? 
is also easily identifiable, and the other two are obtained 
via pseudorotations) . The positive and negative ions are 
slightly triaxial, and as a result their PES's exhibit 12 
local minima, six of them close to oblate shapes and the 
other six close to prolate shapes; in each case, the six 
minima are energetically degenerate and are related to 
each other via pseudorotations. 

At T = 300 K, the canonical-ensemble distribution of 
the initial 10^ phase-space points on the PES of Ag^4 is 
distributed about the 0-like minima, since for the anion 
these minima are the global ones by k, 0.06 eV (or 700 
K) compared to the P-like minima. As a result of the 
vertical photodetachment, some of the initial shapes will 
land on the basins of the P minimum of the neutral Agi4; 
nevertheless most of the initial shapes will land on the 
basins of the O minimum. For the neutral Agi4, these 
latter basins are deep, separated from the basins of the 
P minimum by a barrier of « 0.35 eV (i.e., 4060 K), and 
thus only initial 0-shapes associated with momenta far 
in the tails of the canonical distribution will be able to 
make excursions outside these basins. Consequently, the 
majority of the (constant-energy) trajectories of Agi4 will 
correspond to shape vibrations around the O minimum 
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(see the Lissajous-type trajectory in the middle panel of 
Fig. 1). 

To investigate how these vibrations can be probed by 
a photoionizing laser pulse of a given frequency, we dis- 
play in Fig. 2 the vIP surface [which is the difference 
between the PES of the cation and that of the neutral 
Agi4 (see Fig. 1)]. It is noteworthy that the local min- 
ima of the neutral Agi4 do not correspond to extrema in 
the vIP's. Indeed the O minimum has a vIP of 5.85 eV 
which is intermediate between the extremal values of 6.0 
eV and 5.7 eV. Consequently, a photoionizing laser pulse 
with photon energy 2hu = 5.85 eV is ideal for probing 
the vibrations about the O minimum, since the ionization 
condition vIP < 2'hu) will be satisfied only for certain seg- 
ments of the constant-energy trajectories, and thus the 
measured photoionization yield will exhibit oscillations 
reflecting the shape vibrations. 

An inspection of the probability 7'(i) (which is dis- 
played in Fig. 3) corroborates indeed the above analysis, 
and in addition allows for determination of the period 
of the shape vibrations around the O minimum of Agi4. 
For T = 300 K (lower dashed-dotted curve labelled 1), 
the probability 7^(i) exhibits large amplitude oscillations 
with a period of « 800 fs. Raising the temperature of 
the initial Ag|j~4 ensemble tends to progressively smear 
out the oscillatory structure in 'P(t). This behavior is a 
consequence of the fact that the distribution of the ini- 
tial phase-space points on the PES of the anion is more 
spread away from the region of the 0-like minima (a rel- 
atively larger amount of initial points is distributed over 
the region of the P-like minima) . As a result, a larger, but 
still not dominant, fraction of constant-energy trajecto- 
ries on the PES of Agi4 is restricted to the regions inside 
the basins of the P minima. For T = 600 K (lower dotted 
curve labelled 2), the effect of such trajectories is to re- 
duce the amplitude of the oscillations of the V{t) curve. 
We also note here that for the range of temperatures con- 
sidered (which covers the experimentally expected range) 
no influence on the period of the oscillations is predicted. 

With a somewhat larger probe-laser frequency (e.g., 
2huj = 5.90 eV), the ionization condition vIP < 2hLU is 
satisfied for longer segments of the constant-energy tra- 
jectories, and as a result the troughs in the V{t) curve 
are now less accentuated compared to the case of 2?iuj = 
5.85 eV (see upper curves labelled 3 and 4 in Fig. 3). 

We remark that the experiments suggested in our pa- 
per, and the NeNePo experiments performed in Ref. 
|l3| , involve electron detachment and ionization of clus- 
ters, and they can be performed with photon energies 
which are away from the plasmon excitations of sil- 
ver clusters, which otherwise would have complicated 
the interpretatiop_of such experiments due to resonant 
photoabsorption.a 

In summary, a molecular dynamics in shape space 
method for simulating the dynamics of low-frequency 
collective modes associated with metal-cluster shapes 
was introduced. These collective modes (adiabatic in 
the Born-Oppenheimer sense) are analogous to the well 



known low-energy isoscalar vibrations of atomic nuclei.BI 
The ionic degrees of freedom were described via a dy- 
namically deformable jellium with inertia parameters de- 
rived from an incompressible, irrotational flow, and the 
shell correction method was used to calculate the elec- 
tronic potential energy surface underlying the dynamics 
of the cluster shape. Our finite-temperature simulations 
of Agi4 and its ions, following the negative-to- neutral-to- 
positive scheme, serve to illustrate the MDSS method and 
demonstrate the potential of pump-and-probe ultrashort 
laser pulses as a spectroscopy of cluster shape vibrations. 

We further remark that the MDSS method can be 
generalized to provide a dynamical theory for metal- 
cluster fission (see Refs. [9(a), 24— 26]). In analogy with 
recent dynamical nuclear fission studies,cZI such a theory 
will provide efficient means for calculating fission rates, 
branching ratios between fission channels, the dynamics 
of fission isomers [9(a), 25], as well as the kinetic energy 
distributions of the fission fragments for a broad range 
of cluster sizes. 
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FIG. 1. PES's of Ag^4 (bottom panel), AgM (middle), 
and Agjj (top), a and b denote dimensionless semi-axes of 
the ellipsoidal shapes (the third semi-axis c — 1/ab; see text 
for details). The contour lines correspond to increments of 
0.1 eV in the total energy of the clusters. The Lissajous-type 
figure (heavy solid line), inside one of the basins associated 
with the oblate isomer of Agi4, portrays one of the 3-ps, con- 
stant-energy trajectories at T = 300 K (see text). 

FIG. 2. vIP surface for Agi4. The contour lines corre- 
spond to increments of 0.1 eV. a and b denote dimension- 
less semi-axes of the ellipsoidal shapes (the third semi-ajcis 
c = 1/ab). 

FIG. 3. The time evolution of the probability P(f) (see 
text) for T = 300 K (dashed-dotted lines) and T = 600 K 
(dotted lines). The two-photon energy of the photoionizing 
probe-laser pulse equals 5.85 eV (lower pair of curves labelled 
1 and 2) and 5.90 eV (upper pair of curves labelled 3 and 4). 
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